close all
k = 1e-2;

p = params;

N = 1000;
z = linspace(0, 1, N)';
h = z(2) - z(1);

M = spalloc(N, N, 4 * N);


ch = cosh(k);
sh = sinh(k);
denom = ch * sh - k;


A(1) = -p.Ma / 2 * (k * sh) / denom;
A(2) = 0;
A(3) = p.Ma / 2 * (k*ch - sh) / denom;
A(4) = -A(1) / k;

A = p.delta * A;

w = @(z) (A(1) * z + A(2)) .* cosh(k*z) + (A(3) * z + A(4)) .* sinh(k*z);
c1bz = @(z) -p.beta * (1 - p.beta) * z;

for i = 1:N
    
    if(i == 1)
        M(i,i) = -2 / h^2 - k^2;
        M(i,i+1) = 2 / h^2;
        M(i,N) = -w(z(i)) * c1bz(z(i));
    elseif (i == N - 1)
        M(i,i) = -2 / h^2 - k^2;
        M(i,i-1) = 1 / h^2;
        M(i,i+1) = 1 / h^2 - w(z(i)) * c1bz(z(i));
    elseif (i == N)
        M(i,i) = -2 / h^2 - k^2 -w(z(i)) * c1bz(z(i));
        M(i,i-1) = 2 / h^2;
    else
        M(i,i) = -2 / h^2 - k^2;
        M(i,i-1) = 1 / h^2;
        M(i,i+1) = 1 / h^2;
        M(i,N) = -w(z(i)) * c1bz(z(i));
    end
end


c = 0.01 * ones(size(z));

dt = 1e-3;

for i = 1:round(p.Tmax / dt)
    c = (speye(N) - dt * M) \ c;
    if (mod(i, 10) == 0)
        plot(z, c);
        ylim([0.009, 0.011]);
        pause(0.05);
    end
end